#################################################################
#                       Out Of Sample Model                     #
#################################################################


################################################################################# 
#                     Load Packages                                            #
################################################################################# 


packages = c("tidyverse","zoo","stargazer", "mice",
             "rstan","zoo","bayesplot","pwr", "gridExtra", "loo")
#Warning! This will install R packages in your system if not yet installed. 
#load or install & load all packages 
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

################################################################################# 
#                     Prepare Data & Set Color Scheme                                             #
################################################################################# 
LEDs_data_f <- read.csv(here::here("replication_data_afs.csv"))

color_scheme_set("gray")

############################################################ 
#                       Stan Models WL                    #
############################################################ 

LEDs_data_wl <- LEDs_data_f %>% subset(select = c(-wdldownes,
                                                -wldownes,
                                                -logLER)) %>% drop_na()

########## SUGGESTED WITH CORRUPTION  ##########

col = c( "relative_corr", #1
         "initally",#2
         "targally", #3
         "v2x_libdem", #4
         "v2x_ex_military", #5
         "concap", #6
         "capasst", #7
         "qualrat", #8
         "terrain", #9
         "straterr", #10
         "strat1",  #11
         "strat2", #12
         "strat3",  #13
         "strat4" #14

) 

stan_mirt <- list(
  K =  as.integer(length(col)),
  X = LEDs_data_wl[,col],
  Y = as.numeric(LEDs_data_wl$wl),
  Nd = nrow(LEDs_data_wl)
)

model_f <- stan(file = "simple_logit.stan", data = stan_mirt,
                chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main <- mcmc_areas(model_f,prob = 0.95, pars = c("beta[14]",
                                                      "beta[13]", "beta[12]","beta[11]","beta[10]",
                                                      "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                      "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                      "beta[1]")) +

  scale_y_discrete(labels=c(
    "beta[14]" =  "Strategy 4",
    "beta[13]" =  "Strategy 3",
    "beta[12]" =  "Strategy 2",
    "beta[11]" =  "Strategy 1",
    "beta[10]" =  "Strategy*Terrain",
    "beta[9]" =  "Terrain",
    "beta[8]" = "Troop quality", #8
    "beta[7]" =  "Alliance Capabilities",
    "beta[6]" =  "Capabilities",
    "beta[5]" =  "Military Executive Dependence",
    "beta[4]" =  "Liberal Democracy", 
    "beta[3]" =  "Target", 
    "beta[2]" =  "Initiation",
    "beta[1]" = "Relative Corruption"
  )) 
coef_main

stargazer(summary(model_f, pars = "beta")$summary)

########## SUGGESTED NO CORRUPTION ##########
col_nocorr = c(         "initally",#2
                        "targally", #3
                        "v2x_libdem", #4
                        "v2x_ex_military", #5
                        "concap", #6
                        "capasst", #7
                        "qualrat", #8
                        "terrain", #9
                        "straterr", #10
                        "strat1",  #11
                        "strat2", #12
                        "strat3",  #13
                        "strat4" #14
                        
         
) 

stan_mirt <- list(
  K =  as.integer(length(col_nocorr)),
  X = LEDs_data_wl[,col_nocorr],
  Y = as.numeric(LEDs_data_wl$wl),
  Nd = nrow(LEDs_data_wl)
)

model_f_nocorr <- stan(file = "simple_logit.stan", data = stan_mirt,
                chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_nocorr <- mcmc_areas(model_f_nocorr,prob = 0.95, pars = c("beta[13]", "beta[12]","beta[11]","beta[10]",
                                                      "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                      "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                      "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" =  "Military Executive Dependence",
    "beta[3]" =  "Liberal Democracy", 
    "beta[2]" =  "Target", 
    "beta[1]" =  "Initiation"
  )) 
coef_main_nocorr

stargazer(summary(model_f, pars = "beta")$summary)

########## Reiter and Stam (2002) CORRUPTION ##########

col_RS02_C = c( "relative_corr", #1
         "polini",#2
         "poltarg", #3
         "init", #4
         "concap", #5
         "capasst", #6
         "qualrat", #7
         "terrain", #8
         "straterr", #9
         "strat1",  #10
         "strat2", #11
         "strat3",  #12
         "strat4" #13
) 

stan_mirt <- list(
  K =  as.integer(length(col_RS02_C)),
  X = LEDs_data_wl[,col_RS02_C],
  Y = as.numeric(LEDs_data_wl$wl),
  Nd = nrow(LEDs_data_wl)
)

model_RS02_C <- stan(file = "simple_logit.stan", data = stan_mirt,
                      chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS02_C <- mcmc_areas(model_RS02_C,prob = 0.95, pars = c( "beta[13]", "beta[12]","beta[11]","beta[10]",
                                                                     "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                     "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                     "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[19]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality",
    "beta[6]" = "Alliance Capabilities", #8
    "beta[5]" = "Capabilities",
    "beta[4]" = "Initiation",
    "beta[3]" = "Politics*Target", 
    "beta[2]" = "Politics*Initiation", 
    "beta[1]" =  "Relative Corruption"
  )) 
coef_main_RS02_C


########## Reiter and Stam (2002) NO CORRUPTION ##########

col_RS02_NC = c("polini",#1
         "poltarg", #2
         "init", #3
         "concap", #4
         "capasst", #5
         "qualrat", #6
         "terrain", #7
         "straterr", #8
         "strat1",  #9
         "strat2", #10
         "strat3",  #11
         "strat4" #12
) 


stan_mirt <- list(
  K =  as.integer(length(col_RS02_NC)),
  X = LEDs_data_wl[,col_RS02_NC],
  Y = as.numeric(LEDs_data_wl$wl),
  Nd = nrow(LEDs_data_wl)
)

model_RS02_NC <- stan(file = "simple_logit.stan", data = stan_mirt,
                chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS02_NC <- mcmc_areas(model_RS02_NC,prob = 0.95, pars = c( "beta[12]","beta[11]","beta[10]",
                                                                     "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                     "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                     "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[12]" =  "Strategy 4",
    "beta[11]" =  "Strategy 3",
    "beta[10]" =  "Strategy 2",
    "beta[9]" =  "Strategy 1",
    "beta[8]" =  "Strategy*Terrain",
    "beta[7]" =  "Terrain",
    "beta[6]" = "Troop quality",
    "beta[5]" = "Alliance Capabilities ", #8
    "beta[4]" = "Capabilities",
    "beta[3]" = "Politics*Target",
    "beta[2]" = "Target",
    "beta[1]" = "Politics*Initiation"
  )) 
coef_main_RS02_NC

  
########## Reiter and Stam (2009) CORRUPTION  ##########

col_RS09_C = c("relative_corr", #1
        "pol21targally",#2
        "initally", #3
        "pcini1old", #4
        "pcini2old", #5
        "concap", #6
        "capasst", #7
        "qualrat", #8
        "terrain", #9
        "straterr", #10
        "strat1",  #11
        "strat2", #12
        "strat3",  #13
        "strat4" #14
) 


stan_mirt <- list(
  K =  as.integer(length(col_RS09_C)),
  X = LEDs_data_wl[,col_RS09_C],
  Y = as.numeric(LEDs_data_wl$wl),
  Nd = nrow(LEDs_data_wl)
)

model_RS09_C <- stan(file = "simple_logit.stan", data = stan_mirt,
                      chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS09_C <- mcmc_areas(model_RS09_C,prob = 0.95, pars = c("beta[14]","beta[13]", 
                                                                  "beta[12]","beta[11]","beta[10]",
                                                                    "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                    "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                    "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[14]" =  "Strategy 4",
    "beta[13]" =  "Strategy 3",
    "beta[12]" =  "Strategy 2",
    "beta[11]" =  "Strategy 1",
    "beta[10]" =  "Strategy*Terrain",
    "beta[9]" =  "Terrain",
    "beta[8]" = "Troop quality", #8
    "beta[7]" =  "Alliance Capabilities",
    "beta[6]" =  "Capabilities",
    "beta[5]" = "Politics (second curvilinear term)",
    "beta[4]" = "Politics (first curvilinear term)", 
    "beta[3]" =  "Initiation", 
    "beta[2]" = "Politics*Target", 
    "beta[1]" =  "Relative Corruption"
  )) 
coef_main_RS09_C


########## Reiter and Stam (2009) NOCORRUPTION ##########
col_RS09_NC = c("pol21targally",#1
        "initally", #2
        "pcini1old", #3
        "pcini2old", #4
        "concap", #5
        "capasst", #6
        "qualrat", #7
        "terrain", #8
        "straterr", #9
        "strat1",  #10
        "strat2", #11
        "strat3",  #12
        "strat4" #13
) 

stan_mirt <- list(
  K =  as.integer(length(col_RS09_NC)),
  X = LEDs_data_wl[,col_RS09_NC],
  Y = as.numeric(LEDs_data_wl$wl),
  Nd = nrow(LEDs_data_wl)
)

model_RS09_NC <- stan(file = "simple_logit.stan", data = stan_mirt,
                     chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS09_NC <- mcmc_areas(model_RS09_NC,prob = 0.95, pars = c("beta[13]", 
                                                                    "beta[12]","beta[11]","beta[10]",
                                                                    "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                    "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                    "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" = "Politics (second curvilinear term)",
    "beta[3]" = "Politics (first curvilinear term)", 
    "beta[2]" = "Initiation", 
    "beta[1]" = "Politics*Target"
  )) 

coef_main_RS09_NC

########## Compare ##########

#Main + Corruption
log_lik <- loo::extract_log_lik(model_f)
loo_Main <- loo(log_lik, moment_match = TRUE,)
print(loo_Main)
pkdf<-data.frame(pk=loo_Main$diagnostics$pareto_k,
                 neff=loo_Main$diagnostics$n_eff,
                 n=1:196)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#Main NO Corruption
log_lik_1 <- loo::extract_log_lik(model_f_nocorr)
loo_MainNC <- loo(log_lik_1, moment_match = TRUE)
print(loo_MainNC)
pkdf1<-data.frame(pk=loo_MainNC$diagnostics$pareto_k,
                  neff=loo_MainNC$diagnostics$n_eff,
                  n=1:196)
pkplot1 <- ggplot(pkdf1, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: All Controls") + theme_bw()

#R&S 2002 Corruption
log_lik <- loo::extract_log_lik(model_RS02_C)
loo_RS02C <- loo(log_lik, moment_match = TRUE,)
print(loo_RS02C)
pkdf<-data.frame(pk=loo_RS02C$diagnostics$pareto_k,
                 neff=loo_RS02C$diagnostics$n_eff,
                 n=1:196)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#R&S 2002 No Corruption
log_lik <- loo::extract_log_lik(model_RS02_NC)
loo_RS02NC <- loo(log_lik, moment_match = TRUE,)
print(loo_RS02NC)
pkdf<-data.frame(pk=loo_RS02NC$diagnostics$pareto_k,
                 neff=loo_RS02NC$diagnostics$n_eff,
                 n=1:196)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#R&S 2009  Corruption
log_lik <- loo::extract_log_lik(model_RS09_C)
loo_RS09C <- loo(log_lik, moment_match = TRUE,)
print(loo_RS09C)
pkdf<-data.frame(pk=loo_RS09C$diagnostics$pareto_k,
                 neff=loo_RS09C$diagnostics$n_eff,
                 n=1:196)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#R&S 2009 No Corruption
log_lik <- loo::extract_log_lik(model_RS09_NC)
loo_RS09NC <- loo(log_lik, moment_match = TRUE,)
print(loo_RS09NC)
pkdf<-data.frame(pk=loo_RS09NC$diagnostics$pareto_k,
                 neff=loo_RS09NC$diagnostics$n_eff,
                 n=1:196)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#Main Comp#
comp <- loo_compare(loo_Main, loo_MainNC, loo_RS02C, loo_RS02NC, loo_RS09C, loo_RS09NC )
comp
-6+ c(-1,1)* 1.96*2.9

############################################################ 
#                       Stan Models wldownes               #
############################################################ 
LEDs_data_wldownes <- LEDs_data_f %>% subset(select = c(-wdldownes,
                                                  -wl,
                                                  -logLER)) %>% drop_na()


########## SUGGESTED WITH CORRUPTION  ##########

col = c( "relative_corr", #1
         "initally",#2
         "targally", #3
         "v2x_libdem", #4
         "v2x_ex_military", #5
         "concap", #6
         "capasst", #7
         "qualrat", #8
         "terrain", #9
         "straterr", #10
         "strat1",  #11
         "strat2", #12
         "strat3",  #13
         "strat4" #14
         
) 

stan_mirt <- list(
  K =  as.integer(length(col)),
  X = LEDs_data_wldownes[,col],
  Y = as.numeric(LEDs_data_wldownes$wldownes),
  Nd = nrow(LEDs_data_wldownes)
)

model_Main2 <- stan(file = "simple_logit.stan", data = stan_mirt,
                chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main2 <- mcmc_areas(model_Main2,prob = 0.95, pars = c("beta[14]",
                                                      "beta[13]", "beta[12]","beta[11]","beta[10]",
                                                      "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                      "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                      "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[14]" =  "Strategy 4",
    "beta[13]" =  "Strategy 3",
    "beta[12]" =  "Strategy 2",
    "beta[11]" =  "Strategy 1",
    "beta[10]" =  "Strategy*Terrain",
    "beta[9]" =  "Terrain",
    "beta[8]" = "Troop quality", #8
    "beta[7]" =  "Alliance Capabilities",
    "beta[6]" =  "Capabilities",
    "beta[5]" =  "Military Executive Dependence",
    "beta[4]" =  "Liberal Democracy", 
    "beta[3]" =  "Target", 
    "beta[2]" =  "Initiation",
    "beta[1]" = "Relative Corruption"
  )) 
coef_main2

stargazer(summary(model_f, pars = "beta")$summary)

########## SUGGESTED NO CORRUPTION ##########
col_nocorr = c("initally",#2
               "targally", #3
               "v2x_libdem", #4
               "v2x_ex_military", #5
               "concap", #6
               "capasst", #7
               "qualrat", #8
               "terrain", #9
               "straterr", #10
               "strat1",  #11
               "strat2", #12
               "strat3",  #13
               "strat4" #14
               
) 

stan_mirt <- list(
  K =  as.integer(length(col_nocorr)),
  X = LEDs_data_wldownes[,col_nocorr],
  Y = as.numeric(LEDs_data_wldownes$wldownes),
  Nd = nrow(LEDs_data_wldownes)
)

model_nocorr2 <- stan(file = "simple_logit.stan", data = stan_mirt,
                       chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_nocorr2 <- mcmc_areas(model_nocorr2,prob = 0.95, pars = c("beta[13]", "beta[12]","beta[11]","beta[10]",
                                                                    "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                    "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                    "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" =  "Military Executive Dependence",
    "beta[3]" =  "Liberal Democracy", 
    "beta[2]" =  "Target", 
    "beta[1]" =  "Initiation"
  )) 


########## Downs (2009) CORRUPTION ##########
col_d09C = c("relative_corr", #1
               "pol21",#2
               "initally", #3
               "targally", #4
               "pol21initally", #5
               "pol21targally", #6
               "concap", #7
               "capasst", #8
               "qualrat", #9
               "terrain", #10
               "straterr", #11
               "strat1",  #12
               "strat2", #13
               "strat3",  #14
               "strat4" #15
               
) 

stan_mirt <- list(
  K =  as.integer(length(col_d09C)),
  X = LEDs_data_wldownes[,col_d09C],
  Y = as.numeric(LEDs_data_wldownes$wldownes),
  Nd = nrow(LEDs_data_wldownes)
)

model_d09C <- stan(file = "simple_logit.stan", data = stan_mirt,
                        chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_d09C <- mcmc_areas(model_d09C,prob = 0.95, pars = c("beta[15]", "beta[14]", "beta[13]", 
                                                         "beta[12]","beta[11]","beta[10]",
                                                         "beta[9]","beta[8]", "beta[7]",
                                                         "beta[6]","beta[5]","beta[4]",
                                                         "beta[3]","beta[2]","beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[15]" ="Strategy 4",
    "beta[14]" ="Strategy 3",
    "beta[13]" =  "Strategy 2",
    "beta[12]" =  "Strategy 1",
    "beta[11]" = "Strategy*Terrain" ,
    "beta[10]" =  "Terrain",
    "beta[9]" =  "Troop quality", #8,
    "beta[8]" =  "Alliance Capabilities",
    "beta[7]" = "Capabilities",
    "beta[6]" = "Politics*Target",
    "beta[5]" = "Polity*Initiation",
    "beta[4]" = "Target", 
    "beta[3]" = "Initiation", #3
    "beta[2]" = "Polity",#2
    "beta[1]" =  "Relative Corruption" 
  )) 
coef_d09C

########## Downs (2009) NO CORRUPTION ##########

col_d09NC = c(
             "pol21",#1
             "initally", #2
             "targally", #3
             "pol21initally", #4
             "pol21targally", #5
             "concap", #6
             "capasst", #7
             "qualrat", #8
             "terrain", #9
             "straterr", #110
             "strat1",  #11
             "strat2", #12
             "strat3",  #13
             "strat4" #14
             
) 

stan_mirt <- list(
  K =  as.integer(length(col_d09NC)),
  X = LEDs_data_wldownes[,col_d09NC],
  Y = as.numeric(LEDs_data_wldownes$wldownes),
  Nd = nrow(LEDs_data_wldownes)
)

model_d09NC <- stan(file = "simple_logit.stan", data = stan_mirt,
                   chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_d09NC <- mcmc_areas(model_d09NC, prob = 0.95, pars = c("beta[14]", "beta[13]", 
                                                         "beta[12]","beta[11]","beta[10]",
                                                         "beta[9]","beta[8]", "beta[7]",
                                                         "beta[6]","beta[5]","beta[4]",
                                                         "beta[3]","beta[2]","beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[14]" ="Strategy 4",
    "beta[13]" ="Strategy 3",
    "beta[12]" =  "Strategy 2",
    "beta[11]" =  "Strategy 1",
    "beta[10]" = "Strategy*Terrain" ,
    "beta[9]" =  "Terrain",
    "beta[8]" =  "Troop quality", #8,
    "beta[7]" =  "Alliance Capabilities",
    "beta[6]" = "Capabilities",
    "beta[5]" = "Politics*Target",
    "beta[4]" = "Polity*Initiation",
    "beta[3]" = "Target", 
    "beta[2]" = "Initiation", #3
    "beta[1]" = "Polity"#2
  )) 
coef_d09NC

########## Compare ##########

#Main + Corruption
log_lik2 <- loo::extract_log_lik(model_Main2)
loo_Main2 <- loo(log_lik2, moment_match = TRUE,)
print(loo_Main2)
pkdf<-data.frame(pk=loo_Main2$diagnostics$pareto_k,
                 neff=loo_Main2$diagnostics$n_eff,
                 n=1:232)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#Main NO Corruption
log_lik_3 <- loo::extract_log_lik(model_nocorr2)
loo_MainNC2 <- loo(log_lik_3, moment_match = TRUE)
print(loo_MainNC2)
pkdf1<-data.frame(pk=loo_MainNC2$diagnostics$pareto_k,
                  neff=loo_MainNC2$diagnostics$n_eff,
                  n=1:232)
pkplot1 <- ggplot(pkdf1, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: All Controls") + theme_bw()

#D 2009 Corruption
log_lik <- loo::extract_log_lik(model_d09C)
loo_D09C <- loo(log_lik, moment_match = TRUE,)
print(loo_D09C)
pkdf<-data.frame(pk=loo_D09C$diagnostics$pareto_k,
                 neff=loo_D09C$diagnostics$n_eff,
                 n=1:232)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#D 2009 No Corruption
log_lik <- loo::extract_log_lik(model_d09NC)
loo_D09NC <- loo(log_lik, moment_match = TRUE,)
print(loo_D09NC)
pkdf<-data.frame(pk=loo_D09NC$diagnostics$pareto_k,
                 neff=loo_D09NC$diagnostics$n_eff,
                 n=1:232)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#Main Comp#
comp2 <- loo_compare(loo_Main2, loo_MainNC2,loo_D09C, loo_D09NC )
comp2
-2.8+ c(-1,1)* 1.96*1.2

############################################################ 
#                       Stan Models logLER                 #
############################################################ 

LEDs_data_logLER <- LEDs_data_f %>% subset(select = c(-wdldownes,
                                                        -wl,
                                                        -wldownes)) %>% drop_na()

########## SUGGESTED WITH CORRUPTION  ##########

col = c( "relative_corr", #1
         "initally",#2
         "targally", #3
         "v2x_libdem", #4
         "v2x_ex_military", #5
         "concap", #6
         "capasst", #7
         "qualrat", #8
         "terrain", #9
         "straterr", #10
         "strat1",  #11
         "strat2", #12
         "strat3",  #13
         "strat4" #14
         
) 


stan_mirt <- list(
  K =  as.integer(length(col)),
  X = LEDs_data_logLER[,col],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_Main3 <- stan(file = "simple_linear.stan", data = stan_mirt,
                    chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main3 <- mcmc_areas(model_Main3,prob = 0.95, pars = c("beta[14]",
                                                           "beta[13]", "beta[12]","beta[11]","beta[10]",
                                                           "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                           "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                           "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[14]" =  "Strategy 4",
    "beta[13]" =  "Strategy 3",
    "beta[12]" =  "Strategy 2",
    "beta[11]" =  "Strategy 1",
    "beta[10]" =  "Strategy*Terrain",
    "beta[9]" =  "Terrain",
    "beta[8]" = "Troop quality", #8
    "beta[7]" =  "Alliance Capabilities",
    "beta[6]" =  "Capabilities",
    "beta[5]" =  "Military Executive Dependence",
    "beta[4]" =  "Liberal Democracy", 
    "beta[3]" =  "Target", 
    "beta[2]" =  "Initiation",
    "beta[1]" = "Relative Corruption"
  )) 
coef_main3

########## SUGGESTED NO CORRUPTION ##########
col_nocorr = c("initally",#2
               "targally", #3
               "v2x_libdem", #4
               "v2x_ex_military", #5
               "concap", #6
               "capasst", #7
               "qualrat", #8
               "terrain", #9
               "straterr", #10
               "strat1",  #11
               "strat2", #12
               "strat3",  #13
               "strat4" #14
               
) 


stan_mirt <- list(
  K =  as.integer(length(col_nocorr)),
  X = LEDs_data_logLER[,col_nocorr],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_nocorr3 <- stan(file = "simple_linear.stan", data = stan_mirt,
                      chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_nocorr3 <- mcmc_areas(model_nocorr3,prob = 0.95, pars = c("beta[13]", "beta[12]","beta[11]","beta[10]",
                                                               "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                               "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                               "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" =  "Military Executive Dependence",
    "beta[3]" =  "Liberal Democracy", 
    "beta[2]" =  "Target", 
    "beta[1]" =  "Initiation"
  )) 
coef_nocorr3

########## Reiter and Stam (2002) CORRUPTION ##########

col_RS02_C2 = c( "relative_corr", #1
                "polini",#2
                "poltarg", #3
                "init", #4
                "concap", #5
                "capasst", #6
                "qualrat", #7
                "terrain", #8
                "straterr", #9
                "strat1",  #10
                "strat2", #11
                "strat3",  #12
                "strat4" #13
) 

stan_mirt <- list(
  K =  as.integer(length(col_RS02_C2)),
  X = LEDs_data_logLER[,col_RS02_C2],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_RS02_C2 <- stan(file = "simple_linear.stan", data = stan_mirt,
                     chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS02_C2 <- mcmc_areas(model_RS02_C2,prob = 0.95, pars = c( "beta[13]", "beta[12]","beta[11]","beta[10]",
                                                                   "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                   "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                   "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality",
    "beta[6]" = "Alliance Capabilities", #8
    "beta[5]" = "Capabilities",
    "beta[4]" = "Initiation",
    "beta[3]" = "Politics*Target", 
    "beta[2]" = "Politics*Initiation", 
    "beta[1]" =  "Relative Corruption"
  )) 
coef_main_RS02_C2


########## Reiter and Stam (2002) NO CORRUPTION ##########

col_RS02_NC2 = c("polini",#1
                "poltarg", #3
                "init", #4
                "concap", #5
                "capasst", #6
                "qualrat", #7
                "terrain", #8
                "straterr", #9
                "strat1",  #10
                "strat2", #11
                "strat3",  #12
                "strat4" #13
) 

stan_mirt <- list(
  K =  as.integer(length(col_RS02_NC2)),
  X = LEDs_data_logLER[,col_RS02_NC2],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_RS02_NC2 <- stan(file = "simple_linear.stan", data = stan_mirt,
                      chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS02_NC2 <- mcmc_areas(model_RS02_NC2,prob = 0.95, pars = c(  "beta[12]","beta[11]","beta[10]",
                                                                     "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                     "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                     "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[12]" =  "Strategy 4",
    "beta[11]" =  "Strategy 3",
    "beta[10]" =  "Strategy 2",
    "beta[9]" =  "Strategy 1",
    "beta[8]" =  "Strategy*Terrain",
    "beta[7]" =  "Terrain",
    "beta[6]" = "Troop quality",
    "beta[5]" = "Alliance Capabilities", #8
    "beta[4]" = "Capabilities",
    "beta[3]" = "Initiation",
    "beta[2]" = "Politics*Target", 
    "beta[1]" = "Politics*Initiation"
  )) 
coef_main_RS02_NC2


########## Reiter and Stam (2009) CORRUPTION  ##########
col_RS09_C2 = c("relative_corr", #1
               "pol21targally",#2
               "initally", #3
               "pcini1old", #4
               "pcini2old", #5
               "concap", #6
               "capasst", #7
               "qualrat", #8
               "terrain", #9
               "straterr", #10
               "strat1",  #11
               "strat2", #12
               "strat3",  #13
               "strat4" #14
) 

stan_mirt <- list(
  K =  as.integer(length(col_RS09_C2)),
  X = LEDs_data_logLER[,col_RS09_C2],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_RS09_C2 <- stan(file = "simple_linear.stan", data = stan_mirt,
                     chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS09_C2 <- mcmc_areas(model_RS09_C2,prob = 0.95, pars = c("beta[14]","beta[13]", 
                                                                  "beta[12]","beta[11]","beta[10]",
                                                                  "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                  "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                  "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[14]" =  "Strategy 4",
    "beta[13]" =  "Strategy 3",
    "beta[12]" =  "Strategy 2",
    "beta[11]" =  "Strategy 1",
    "beta[10]" =  "Strategy*Terrain",
    "beta[9]" =  "Terrain",
    "beta[8]" = "Troop quality", #8
    "beta[7]" =  "Alliance Capabilities",
    "beta[6]" =  "Capabilities",
    "beta[5]" = "Politics (second curvilinear term)",
    "beta[4]" = "Politics (first curvilinear term)", 
    "beta[3]" = "Initiation", 
    "beta[2]" = "Politics*Target", 
    "beta[1]" =  "Relative Corruption"
  )) 
coef_main_RS09_C2


########## Reiter and Stam (2009) NOCORRUPTION ##########
col_RS09_NC2 = c("pol21targally",#1
                "initally", #2
                "pcini1old", #3
                "pcini2old", #4
                "concap", #5
                "capasst", #6
                "qualrat", #7
                "terrain", #8
                "straterr", #9
                "strat1",  #10
                "strat2", #11
                "strat3",  #12
                "strat4" #13
) 

stan_mirt <- list(
  K =  as.integer(length(col_RS09_NC2)),
  X = LEDs_data_logLER[,col_RS09_NC2],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_RS09_NC2 <- stan(file = "simple_linear.stan", data = stan_mirt,
                      chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_RS09_NC2 <- mcmc_areas(model_RS09_NC2,prob = 0.95, pars = c("beta[13]", 
                                                                    "beta[12]","beta[11]","beta[10]",
                                                                    "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                    "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                    "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" = "Politics (second curvilinear term)",
    "beta[3]" = "Politics (first curvilinear term)", 
    "beta[2]" = "Initiation", 
    "beta[1]" = "Politics*Target"
  )) 

coef_main_RS09_NC2


########## Downs (2009) CORRUPTION ##########
col_d09C = c("relative_corr", #1
             "pol21",#2
             "initally", #3
             "targally", #4
             "pol21initally", #5
             "pol21targally", #6
             "concap", #7
             "capasst", #8
             "qualrat", #9
             "terrain", #10
             "straterr", #11
             "strat1",  #12
             "strat2", #13
             "strat3",  #14
             "strat4" #15
             
) 

stan_mirt <- list(
  K =  as.integer(length(col_d09C)),
  X = LEDs_data_logLER[,col_d09C],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_d09C2 <- stan(file = "simple_linear.stan", data = stan_mirt,
                   chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_d09C2 <- mcmc_areas(model_d09C2,prob = 0.95, pars = c("beta[15]", "beta[14]", "beta[13]", 
                                                         "beta[12]","beta[11]","beta[10]",
                                                         "beta[9]","beta[8]", "beta[7]",
                                                         "beta[6]","beta[5]","beta[4]",
                                                         "beta[3]","beta[2]","beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[15]" ="Strategy 4",
    "beta[14]" ="Strategy 3",
    "beta[13]" =  "Strategy 2",
    "beta[12]" =  "Strategy 1",
    "beta[11]" = "Strategy*Terrain" ,
    "beta[10]" =  "Terrain",
    "beta[9]" =  "Troop quality", #8,
    "beta[8]" =  "Alliance Capabilities",
    "beta[7]" = "Capabilities",
    "beta[6]" = "Politics*Target",
    "beta[5]" = "Polity*Initiation",
    "beta[4]" = "Target", 
    "beta[3]" = "Initiation", #3
    "beta[2]" = "Polity",#2
    "beta[1]" =  "Relative Corruption" 
  )) 
coef_d09C2

########## Downs (2009) NO CORRUPTION ##########

col_d09NC = c(
  "pol21",#1
  "initally", #2
  "targally", #3
  "pol21initally", #4
  "pol21targally", #5
  "concap", #6
  "capasst", #7
  "qualrat", #8
  "terrain", #9
  "straterr", #110
  "strat1",  #11
  "strat2", #12
  "strat3",  #13
  "strat4" #14
  
) 

stan_mirt <- list(
  K =  as.integer(length(col_d09NC)),
  X = LEDs_data_logLER[,col_d09NC],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_d09NC2 <- stan(file = "simple_linear.stan", data = stan_mirt,
                    chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_d09NC2 <- mcmc_areas(model_d09NC2, prob = 0.95, pars = c("beta[14]", "beta[13]", 
                                                            "beta[12]","beta[11]","beta[10]",
                                                            "beta[9]","beta[8]", "beta[7]",
                                                            "beta[6]","beta[5]","beta[4]",
                                                            "beta[3]","beta[2]","beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[14]" ="Strategy 4",
    "beta[13]" ="Strategy 3",
    "beta[12]" =  "Strategy 2",
    "beta[11]" =  "Strategy 1",
    "beta[10]" = "Strategy*Terrain" ,
    "beta[9]" =  "Terrain",
    "beta[8]" =  "Troop quality", #8,
    "beta[7]" =  "Alliance Capabilities",
    "beta[6]" = "Capabilities",
    "beta[5]" = "Politics*Target",
    "beta[4]" = "Polity*Initiation",
    "beta[3]" = "Target", 
    "beta[2]" = "Initiation", #3
    "beta[1]" = "Polity"#2
  )) 
coef_d09NC2

########## Compare ##########

#Main + Corruption
log_lik5 <- loo::extract_log_lik(model_Main3)
loo_Main3 <- loo(log_lik5, moment_match = TRUE,)
print(loo_Main3)
pkdf<-data.frame(pk=loo_Main3$diagnostics$pareto_k,
                 neff=loo_Main3$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#Main NO Corruption
log_lik_6 <- loo::extract_log_lik(model_nocorr3)
loo_MainNC3 <- loo(log_lik_6, moment_match = TRUE)
print(loo_MainNC3)
pkdf1<-data.frame(pk=loo_MainNC3$diagnostics$pareto_k,
                  neff=loo_MainNC3$diagnostics$n_eff,
                  n=1:192)
pkplot1 <- ggplot(pkdf1, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: All Controls") + theme_bw()

#R&S 2002 Corruption
log_lik <- loo::extract_log_lik(model_RS02_C2)
loo_RS02C2 <- loo(log_lik, moment_match = TRUE,)
print(loo_RS02C2)
pkdf<-data.frame(pk=loo_RS02C2$diagnostics$pareto_k,
                 neff=loo_RS02C2$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#R&S 2002 No Corruption
log_lik <- loo::extract_log_lik(model_RS02_NC2)
loo_RS02NC2 <- loo(log_lik, moment_match = TRUE,)
print(loo_RS02NC2)
pkdf<-data.frame(pk=loo_RS02NC2$diagnostics$pareto_k,
                 neff=loo_RS02NC2$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#R&S 2009  Corruption
log_lik <- loo::extract_log_lik(model_RS09_C2)
loo_RS09C2 <- loo(log_lik, moment_match = TRUE,)
print(loo_RS09C2)
pkdf<-data.frame(pk=loo_RS09C2$diagnostics$pareto_k,
                 neff=loo_RS09C2$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#R&S 2009 No Corruption
log_lik <- loo::extract_log_lik(model_RS09_NC2)
loo_RS09NC2 <- loo(log_lik, moment_match = TRUE,)
print(loo_RS09NC2)
pkdf<-data.frame(pk=loo_RS09NC2$diagnostics$pareto_k,
                 neff=loo_RS09NC2$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()


#D 2009 Corruption
log_lik <- loo::extract_log_lik(model_d09C2)
loo_D09C2 <- loo(log_lik, moment_match = TRUE,)
print(loo_D09C2)
pkdf<-data.frame(pk=loo_D09C2$diagnostics$pareto_k,
                 neff=loo_D09C2$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#D 2009 No Corruption
log_lik <- loo::extract_log_lik(model_d09NC2)
loo_D09NC2 <- loo(log_lik, moment_match = TRUE,)
print(loo_D09NC2)
pkdf<-data.frame(pk=loo_D09NC2$diagnostics$pareto_k,
                 neff=loo_D09NC2$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()




#Main Comp#
comp3 <- loo_compare(loo_Main3, loo_MainNC3,loo_RS02C2, loo_RS02NC2, loo_RS09C2, loo_RS09NC2, loo_D09C2, loo_D09NC2)
comp3

######################################################
######### Print Comparisons for Table 2 #########
######################################################

#RS2002 WL
stargazer(comp)
#calculate if diff is significant to model without corruption
loo_compare(loo_Main, loo_MainNC, loo_RS02C, loo_RS02NC, loo_RS09C, loo_RS09NC )
#prop
print(loo_compare(loo_Main, loo_MainNC), simplify = F)
#RS2002
print(loo_compare(loo_RS02C, loo_RS02NC), simplify = F)
#RS2009
print(loo_compare(loo_RS09C, loo_RS09NC), simplify = F)

#D2009WL
stargazer(comp2)
#calculate if diff is significant to model without corruption
loo_compare(loo_Main2, loo_MainNC2,loo_D09C, loo_D09NC )
#prop
print(loo_compare(loo_Main2, loo_MainNC2), simplify = F)
#D2009
print(loo_compare(loo_D09C, loo_D09NC ), simplify = F)

#LogLER
stargazer(comp3)
loo_compare(loo_Main3, loo_MainNC3,loo_RS02C2, loo_RS02NC2, loo_RS09C2, loo_RS09NC2, loo_D09C2, loo_D09NC2)
#calculate if diff is significant to model without corruption
#prop
print(loo_compare(loo_Main3, loo_MainNC3), simplify = F)
#RS2002
print(loo_compare(loo_RS02C2, loo_RS02NC2), simplify = F)
#RS2009
print(loo_compare(loo_RS09C2, loo_RS09NC2), simplify = F)
#D2009
print(loo_compare(loo_D09C2, loo_D09NC2), simplify = F)     


#Print Tables for Appendix 
###Combine Plots
coef_main<- coef_main + labs(title = "Corruption Model")
coef_main_nocorr<- coef_main_nocorr + labs(title = "No Corruption Model")
grid.arrange(coef_main,coef_main_nocorr,
             layout_matrix = rbind(c(1,2)
             ))

coef_main_RS02_C<- coef_main_RS02_C+ labs(title = "Corruption Model")
coef_main_RS02_NC<- coef_main_RS02_NC+ labs(title = "No Corruption Model")
grid.arrange(coef_main_RS02_C,
             coef_main_RS02_NC,
             layout_matrix = rbind(c(1,2)
             ))

coef_main_RS09_C<- coef_main_RS09_C + labs(title = "Corruption Model")
coef_main_RS09_NC<- coef_main_RS09_NC+ labs(title = "No Corruption Model")
grid.arrange(coef_main_RS09_C,
             coef_main_RS09_NC,
             layout_matrix = rbind(c(1,2)
             ))

coef_main2<- coef_main2+ labs(title = "Corruption Model")
coef_nocorr2<- coef_nocorr2+ labs(title = "No Corruption Model")
grid.arrange(coef_main2,
             coef_nocorr2,
             layout_matrix = rbind(c(1,2)
             ))

coef_d09C<- coef_d09C+ labs(title = "Corruption Model")
coef_d09NC<- coef_d09NC+ labs(title = "No Corruption Model")
grid.arrange(coef_d09C,
             coef_d09NC,
             layout_matrix = rbind(c(1,2)
             ))


coef_main3<- coef_main3 + labs(title = "Corruption Model")
coef_nocorr3<- coef_nocorr3+ labs(title = "No Corruption Model")
grid.arrange(coef_main3,
             coef_nocorr3,
             layout_matrix = rbind(c(1,2)
             ))

coef_main_RS02_C2<- coef_main_RS02_C2+ labs(title = "Corruption Model")
coef_main_RS02_NC2<- coef_main_RS02_NC2+ labs(title = "No Corruption Model")
grid.arrange(coef_main_RS02_C2,
             coef_main_RS02_NC2,
             layout_matrix = rbind(c(1,2)
             ))

coef_main_RS09_C2<- coef_main_RS09_C2+ labs(title = "Corruption Model")
coef_main_RS09_NC2<- coef_main_RS09_NC2+ labs(title = "No Corruption Model")
grid.arrange(coef_main_RS09_C2,
             coef_main_RS09_NC2,
             layout_matrix = rbind(c(1,2)
             ))

coef_d09C2<- coef_d09C2 + labs(title = "Corruption Model")
coef_d09NC2<- coef_d09NC2 + labs(title = "No Corruption Model")        
grid.arrange(coef_d09C2,
             coef_d09NC2,
             layout_matrix = rbind(c(1,2)
             ))

########## SUGGESTED WITH CORRUPTION BUT NO LIBDEM  ##########
#check how much extra information does liberal democracy add

###########Stan Models WL##########

########## SUGGESTED WITH CORRUPTION  ##########

col = c( "relative_corr", #1
         "initally",#2
         "targally", #3
         "v2x_ex_military", #5
         "concap", #6
         "capasst", #7
         "qualrat", #8
         "terrain", #9
         "straterr", #10
         "strat1",  #11
         "strat2", #12
         "strat3",  #13
         "strat4" #14
         
) 

stan_mirt <- list(
  K =  as.integer(length(col)),
  X = LEDs_data_wl[,col],
  Y = as.numeric(LEDs_data_wl$wl),
  Nd = nrow(LEDs_data_wl)
)

model_f_nld <- stan(file = "simple_logit.stan", data = stan_mirt,
                    chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main_nld <- mcmc_areas(model_f_nld,prob = 0.95, pars = c("beta[13]", "beta[12]","beta[11]","beta[10]",
                                                              "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                              "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                              "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" =  "Military Executive Dependence",
    "beta[3]" =  "Target", 
    "beta[2]" =  "Initiation",
    "beta[1]" = "Relative Corruption"
  )) 
coef_main_nld

stargazer(summary(model_f_nld, pars = "beta")$summary)

##########Stan Models wldownes########

col = c( "relative_corr", #1
         "initally",#2
         "targally", #3
         "v2x_ex_military", #4
         "concap", #5
         "capasst", #6
         "qualrat", #7
         "terrain", #8
         "straterr", #9
         "strat1",  #10
         "strat2", #11
         "strat3",  #12
         "strat4" #13
         
) 

stan_mirt <- list(
  K =  as.integer(length(col)),
  X = LEDs_data_wldownes[,col],
  Y = as.numeric(LEDs_data_wldownes$wldownes),
  Nd = nrow(LEDs_data_wldownes)
)

model_Main2_nld <- stan(file = "simple_logit.stan", data = stan_mirt,
                        chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main2_nld <- mcmc_areas(model_Main2_nld,prob = 0.95, pars = c("beta[13]", "beta[12]","beta[11]","beta[10]",
                                                                   "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                   "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                   "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" =  "Military Executive Dependence",
    "beta[3]" =  "Target", 
    "beta[2]" =  "Initiation",
    "beta[1]" = "Relative Corruption"
  )) 
coef_main2_nld

stargazer(summary(model_f_nld, pars = "beta")$summary)


##########Stan Models logLER  ##########

col = c( "relative_corr", #1
         "initally",#2
         "targally", #3
         "v2x_ex_military", #5
         "concap", #6
         "capasst", #7
         "qualrat", #8
         "terrain", #9
         "straterr", #10
         "strat1",  #11
         "strat2", #12
         "strat3",  #13
         "strat4" #14
         
) 


stan_mirt <- list(
  K =  as.integer(length(col)),
  X = LEDs_data_logLER[,col],
  Y = as.numeric(LEDs_data_logLER$logLER),
  Nd = nrow(LEDs_data_logLER)
)

model_Main3_nld <- stan(file = "simple_linear.stan", data = stan_mirt,
                        chains = 2, iter =3000, cores = 2, seed = 2023, control = list(max_treedepth = 15, adapt_delta = 0.95))


coef_main3_nld <- mcmc_areas(model_Main3_nld,prob = 0.95, pars = c("beta[13]", "beta[12]","beta[11]","beta[10]",
                                                                   "beta[9]","beta[8]", "beta[7]","beta[6]",
                                                                   "beta[5]","beta[4]", "beta[3]","beta[2]",
                                                                   "beta[1]")) +
  
  scale_y_discrete(labels=c(
    "beta[13]" =  "Strategy 4",
    "beta[12]" =  "Strategy 3",
    "beta[11]" =  "Strategy 2",
    "beta[10]" =  "Strategy 1",
    "beta[9]" =  "Strategy*Terrain",
    "beta[8]" =  "Terrain",
    "beta[7]" = "Troop quality", #8
    "beta[6]" =  "Alliance Capabilities",
    "beta[5]" =  "Capabilities",
    "beta[4]" =  "Military Executive Dependence",
    "beta[3]" =  "Target", 
    "beta[2]" =  "Initiation",
    "beta[1]" = "Relative Corruption"
  )) 
coef_main3_nld

#Print Tables  
###Combine Plots
coef_main <- coef_main + labs(title = "Liberal Democracy + Corruption ")
coef_main_nocorr <- coef_main_nocorr + labs(title = "Liberal Democracy ")
coef_main_nld <- coef_main_nld + labs(title = "Corruption")
grid.arrange(grobs = list(coef_main,coef_main_nocorr,coef_main_nld),
             layout_matrix = rbind(c(1,2,3)
             ))

coef_main2 <- coef_main2 + labs(title = "Liberal Democracy + Corruption ")
coef_nocorr2 <- coef_nocorr2 + labs(title = "Liberal Democracy ")
coef_main2_nld <- coef_main2_nld + labs(title = "Corruption")
grid.arrange(grobs = list(coef_main2,coef_nocorr2,coef_main2_nld),
             layout_matrix = rbind(c(1,2,3)
             ))

coef_main3 <- coef_main2 + labs(title = "Liberal Democracy + Corruption ")
coef_nocorr3 <- coef_nocorr3 + labs(title = "Liberal Democracy ")
coef_main3_nld <- coef_main3_nld + labs(title = "Corruption")
grid.arrange(grobs = list(coef_main3,coef_nocorr3,coef_main3_nld),
             layout_matrix = rbind(c(1,2,3)
             ))

########Compare#######
#WLR&S
log_lik_f_nld <- loo::extract_log_lik(model_f_nld)
loo_f_nld <- loo(log_lik_f_nld, moment_match = TRUE,)
print(loo_f_nld)
pkdf<-data.frame(pk=loo_f_nld$diagnostics$pareto_k,
                 neff=loo_f_nld$diagnostics$n_eff,
                 n=1:196)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#WLD09
log_lik_Main2_nld <- loo::extract_log_lik(model_Main2_nld)
loo_Main2_nld <- loo(log_lik_Main2_nld, moment_match = TRUE,)
print(loo_Main2_nld)
pkdf<-data.frame(pk=loo_Main2_nld$diagnostics$pareto_k,
                 neff=loo_Main2_nld$diagnostics$n_eff,
                 n=1:232)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

#LogLer
log_lik_Main3_nld <- loo::extract_log_lik(model_Main3_nld)
loo_Main3_nld <- loo(log_lik_Main3_nld, moment_match = TRUE,)
print(loo_Main3_nld)
pkdf<-data.frame(pk=loo_Main3_nld$diagnostics$pareto_k,
                 neff=loo_Main3_nld$diagnostics$n_eff,
                 n=1:192)
pkplot <- ggplot(pkdf, aes(x=n,y=pk)) + geom_point(shape=3, color="blue" ,size=1.5) +
  labs(x="Observation left out", y="Pareto k") +
  geom_hline(yintercept = 0, linetype=3, color="red", size=0.2) +
  geom_hline(yintercept = 0.5, linetype=4, color="red", size=0.2) +
  geom_hline(yintercept = 0.7, linetype=2, color="red", size=0.2) +
  scale_y_continuous(breaks=c(-.25,0, 0.25, 0.5, 0.7,1), limits=c(-0.3,1)) +
  ggtitle("PSIS-LOO diagnostics: Full Model") + theme_bw()

####Tables#####

#WLR&S
print(loo_compare(loo_Main, loo_MainNC, loo_f_nld),simplify = F)
print(loo_compare(loo_Main, loo_f_nld),simplify = F)
print(loo_compare(loo_Main, loo_MainNC),simplify = F)
print(loo_compare(loo_MainNC, loo_f_nld),simplify = F)

#WLD
print(loo_compare(loo_Main2, loo_MainNC2, loo_Main2_nld),simplify = F)
print(loo_compare(loo_Main2, loo_Main2_nld),simplify = F)
print(loo_compare(loo_Main2, loo_MainNC2),simplify = F)
print(loo_compare(loo_MainNC2, loo_Main2_nld),simplify = F)

#LogLer
print(loo_compare(loo_Main3, loo_MainNC3, loo_Main3_nld),simplify = F)
print(loo_compare(loo_Main3, loo_Main3_nld),simplify = F)
print(loo_compare(loo_Main3, loo_MainNC3),simplify = F)
print(loo_compare(loo_MainNC3, loo_Main3_nld),simplify = F)
